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Abstract 

We study numerical methods for the solution of general linear moment prob- 
lems, where the solution belongs to a family of nested subspaces of a Hilbert 
space. Multi-level algorithms, based on the conjugate gradient method and the 
Landweber-Richardson method are proposed that determine the "optimal" recon- 
struction level a posteriori from quantities that arise during the numerical calcu- 
lations. As an important example we discuss the reconstruction of band-limited 
signals from irregularly spaced noisy samples, when the actual bandwidth of the 
signal is not available. Numerical examples show the usefulness of the proposed 
algorithms. 
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1 Introduction 

We study numerical methods for the solution of a family of general linear moment prob- 
lems 

(l.l) (x,g?)=vi j jel 



where {gf}j^i is a complete system in X^ and {X^jei is a family of nested subspaces 
of a separable Hilbert space X. (x,g^) denotes the j—th moment. We consider the 
problem of recovering x from noisy observations of the moments fij + 5j. 

Most existing methods for moments problems require exact a priori information about 
the space Xn to which x belongs to | BG68| , [XN94 , LM91 |. In many applications however 



this knowledge is not available, often we only know that x belongs to some space of a 
whole family of Hilbert spaces X^. 

An important example is the problem of reconstructing a band-limited signal x from 
irregularly spaced, noisy samples x(tj) + Sj, without explicit knowledge on the bandwidth 
of x. In other words we know that x belongs to a space of a family of band-limited 
functions X^ (where N represents the bandwidth), but we do not know to which one. 
Most reconstruction algorithms require the a priori information of the actual bandwidth 
of x (or at least a very good estimate) to be applicable. In numerical reconstructions 
overestimating the bandwidth of x may result in a highly oscillating function, due to the 
noise in the samples. Underestimating the bandwidth may give a poor approximation 
to the original signal. From this point of view a proper choice of the level N can also be 
seen as a regularization procedure. 

For the numerical solution of the family of general moment problems we propose a multi- 
level approach that determines an "optimal" level N a-posteriori from quantities that 
arise during the numerical calculations. The proposed algorithm does not use a-priori 
information on the level N of the solution - all information required in the analysis of 
this algorithm is that there exists a level (bandwidth) N* such that the solution x of 
( |1 . 1| ) satisfies G X^. 

The outline of this paper is the following: In Section || and Section |3| we study multi- 
level algorithms for the efficient stable solution of moment problems. In particular we 
concentrate on a multi- level conjugate gradient algorithm and a multi- level Landweber- 
Richardson algorithm. These algorithms are formulated in a general setting. In Section 
H several applications to the numerical solution of moment problems are given, with the 
emphasize on reproducing kernel Hilbert spaces. In Section |5| we demonstrate in detail 
how the proposed algorithms can be applied to the problem of reconstructing a band- 
limited signal of unknown bandwidth from noisy irregularly spaced samples. Finally 
in Section |6] we demonstrate the performance of the multi-level algorithms by some 
numerical examples. 
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2 A multi— level algorithm based on the conjugate 
gradient method 

In this section we introduce a multi-level algorithm based on the conjugate gradient for 
the numerical solution of moment problems. In this algorithm the levels are adapted 
inductively and at each level a conjugate gradient method is implemented. 

In order to make the paper self-contained we outline the conjugate gradient method for 
the solution of linear systems: Let T : X — > Y be a linear operator between separable 
Hilbert spaces X and Y. By T* we denote the adjoint of the operator T. One efficient 
method for the solution of Tx = y is the conjugate gradient method applied to the normal 
equation T*xT = T*y. This method, which is denoted by CGNE in the literature, is 
defined as follows (see e.g. [ |Han95| , [Han96| , [fcCHJN96|| ) 
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Before we study a multi-level scheme based on the CGNE method we introduce some 
notation: 

Assumption 2.1 {^C/v}jv 6 pj a family of linearly ordered subspaces of a separable 
Hilbert space X . For N £ No let Pn be the orthogonal projectors from X into Xn 
and let TZ(T^) denote the range of T^. For N £ Nq the family of linear operators 
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T/v : Xn — > K satisfies 



W M<1, 



(2.1) <| (ii) T N x M = T M x M forallN>M, 

(iii) ^(T^)CX,v. 

The first assumption in ( [2.1| ) is that the operators Tjy are properly scaled. The second 
assumption is that on subspaces of level M(< iV) the operators Tm and Tjy coincide. 
The third assumptions guarantees that the iterates of the CGNE method are contained 
in X N . 

For given data y £ Y we are concerned with finding JVeN and x G Xjy, which satisfy 

(2.2) T N x = y . 



The terminology of a solution of ( |2.2| ) is well-defined by Assumption ( |2.1| ), since from 
this it follows that a solution on a lower scale is also a solution on a higher scale. 

In the sequel we will assume that there exists a solution of (|2.2|) . 

Assumption 2.2 There exists A* e No and a;* e Ajv„ which satisfy $TQ). 

We note that we do not assume the knowledge of A*. If for instance X^v denotes the 
Payley- Wiener space of band-limited functions of bandwidth less or equal A and A — > 



oo, then the meaning of Assumption [2.2| is: There exists a band-limited function which 
solves (O). 



By ^ 6 F we denote measured data of ?/, which satisfies 

(2-3) \\y s -y\\<s. 



Algorithm 2.3 Let rj > 0, r > 1. 

/eve/ : Without loss of generality we assume that Xq E Xq satisfies 
1. \\T (x 5 0fi )-y 5 \\ 2 > 2(1 + 7?) (5+\\P x*-x m \\) \\y s \\, 



2. ||T « )- 2 / <5 || 2 >2(l + r / )r5||^ 
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does not hold, then Xq is accepted as an approximate solution and the multi- 
level algorithm is terminated. If^. does not hold, but ^. is satisfied, then we can 
choose a level J such that 

\\Tj{x%) - y 5 f = ||To« ) - /f > 2(1 + rj) (6 + \\Pjx* - x m \\) \\y s \\ . 

In this situation we simply renumber the levels and make the setting = J. For 
the new setting J = both [Z|. and |^. hold, and the CGNE method at level is 
performed. The iteration at scale X is terminated if for the first time 

\\To(4+i,o) - V S f < 2(1 + V) (S + \\Pox. -x.\\) \\V 5 \\ ■ 
The termination index is denoted by k = fc*(0, 8). 

From level N to level N + 1: If \\T N+1 (x 5 h ^ NjS)+ltN ) - y 5 \\ 2 < 2(1 + r])T5\\y 5 \\, then the 
multi-level algorithm terminates and a^ny^+i n ^ s an approximate solution. Else 
we define x s k , N ^ N =: Xq >n+1 . Since HPjv+i^* — ^*|| < ll-P/v^* — ^*|| it follows from 
the definition of k*(N,6): 

\\T N+1 {x s 0tN+1 ) - y s \\ > 2(1 + rj) (5 + \\P N+1 x* - x4) \\y s \\ . 

Therefore, at least one iteration of CGNE at level N+l is performed. The iteration 
at scale X^+i is terminated if 

\\Tn+i(xI + i,n+i) ~ y S \\ 2 < 2(1 + V) (8 + \\P N+ ix* -x # \\) \\y s \\ 
for the first time. The termination index is denoted by k = k*(N + 1,5). 

The idea of using multi-level algorithms as considered in this paper for the solution 
of ill-posed problems has been considered in | jjch96|| . 

In the following, we verify some basic monotonicity properties of the iterates of 
the multi-level algorithm |2.3| . Note that in practice when the CGNE-method is used 
for solving linear systems the CGNE method is usually terminated with a discrepancy 



principle (see e.g. |Han96|). Hanke [ Han96 | showed by a counterexample that in general 



the error of the iterates of the CGNE method does not monotonically decrease when the 
iteration is terminated by the discrepancy principle. In this paper the CGNE-method 
is combined with a multi-level algorithm and we expect that a reasonable performance 
of this method is depending on providing a good initially guess for starting the CGNE 
method at a higher level. Therefore, we concentrate in this paper on a stopping criterion, 
which guarantees monotonicity of the iterates at each level. 
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Proposition 2.4 Let assumption \2.1\ hold, and let the iteration at fixed level N be 
stopped according to the generalized discrepancy principle 

(2.4) \\T N (xi >N ) - y 5 \\ 2 > 2(1 + 77) (<S + \\P N x. - \\y 5 \\ forO<k< K(N, 5) . 
Then (N, 5) < 00 and 

\\xf +1 - P N x*\\ < \\xf - Pnx*\\ j = 0, h(N, 8) , 



- 1 1 -/V * 1 1 _ \\ 2 

xf +1 — a;* 1 1 < H^-x*!! j = 0, K(N, 5) . 



Proof: Application of the general results in ||Han96|| immediately yields (N, 5) < 00 



and for j = 0, k*(N, 5) — %*\\ < — Therefore we have 

11 JV _ p#~ ||2 , 9 /,JV _ p N p N _ \ < II N_ „ ||2 , / JV _ pN pN T _ T \ 

+ 1 ** / * 1 1 \ J ' I 1 * j ^ * ** / * / _ 1 1 J * 1 1 \ j I 1 * ? * 



From the remark following Assumption [O] it follows Xj ,Xj +x G Xjv and therefore the 



assertion is proven. □ 



From Proposition |2]4] it follows that fl2.4|) determines a well-defined stopping index for 



the CGNE method at level N. 

2.1 Convergence and stability of the multi— level algorithm 

Using the monotonicity results of the iterates of the multi-level schemes we prove below 
that the multi-level scheme is convergent and stable. The following lemma will be central 
for our further considerations. 



Lemma 2.5 Let 5 = 0. If Assumptions 2.1 and 2.& holds, then the multi-level algorithm 



2.5 



terminates after a finite number of iterations 



or 



there exists a level N such that at each level N < N a finite number of CGNE-iterations 
are performed and at level N the CGNE iteration does not terminate. 



Proof: If 5 = 0, then the stopping criterion (O) for the CGNE method at level N is 

(2.5) \\T n (4,n) ~yf> 2(1 + V) \\Pnx* ~ x m \\ \\y\\ for < k < K{N, 0) . 
If 5 = 0, then the stopping criterion for the multi-level algorithm is 

(2.6) \\T N {xi N )-y\\ = Q. 
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Let N be the first level-index where A finite level-index exists by As- 

sumption |2.2|. From Proposition [2.4| it follows that for each level N < N the number of 



iterations of the CGNE method is finite, which yields the assertion. □ 

Using the previous lemma, we are able to prove that the multi-level CGNE algorithm is 
convergent in the case of exact data. 

Theorem 2.6 Let Assumption \2.2\ hold. Then the multi-level CGNE algorithm con- 
verges to a solution x. 

Proof: Let iV be the smallest index N such that P^x* = x*. We differ between two 
cases: 

1. The multi-level CGNE method terminates at level N after a finite number of 
iterations 

2. The multi-level CGNE method does not terminate 

In the first case it follows from Q2.6| ) that the last iterate x satisfies T^x = y , which 
shows that x is the solution. 



In the second case it follows immediately from Theorem 3.4 in |Han95] that x n is con- 



vergent to a solution. □ 
In the case of perturbed data, it follows from Proposition |2.4| that the discrepancy 



principle ( |2.4| ) determines a well-defined (finite) stopping index at level N < N, which 
is denoted by k*(N,8). In the following, we will assume that the multi-level algorithm 
is finally terminated if for the first time 



(2-7) \\T N (x{ N )-y s f < 2(1 + V )T8\\y 



From Lemma 2.5 and Assumption [2.2| it follows that this criterion determines a well- 



defined stopping index at level N(S), which is denoted by K*(N(6),6). 
Our next result shows that the multi-level iteration is a stable method. 



Theorem 2.7 Let the assumptions of Theorem \2. (\ hold and assume that Tjv. is contin- 
uously invertible (here N* is as introduced in Assumption \2.<\ ). If y 5 fulfills ( \2. 3j) , and if 
the iteration at level N is stopped according to the discrepancy principle ($-4\ ), and the 
multi-level iteration is terminated by \2. then 



J - 5^0 . 



X N n ,K*{N n ,8) ~^ x *i 
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Proof: Let 8 n , n = 1,2,..., be a sequence converging to zero as n — > oo, and let 



Vn '■= y 5n be a corresponding sequence of perturbed data. For each pair (5 n , y n ), denote by 
K n = K^NlS^jSn) the corresponding stopping index at termination level N n = N(S n ) 
determined by (|2.?1 ). A comparison of 5 n + \\Pn(s„)X* — £*|| and r5 n shows that N(5 n ) < 
iV*. From Proposition O and the results in [|Han9(: | it follows that for < k < k*(N, S) 
and < N < N* 



\ X k,N X * 



< 



\ x X * I 



and \\T Nt x 5 kN - y\ 







From the continuous invertibility of the assertion follows. 



□ 



3 A multi— level algorithm based on the Landweber— 
Richardson method 

As an alternative to the multi-level algorithm based on CGNE we study in this section 
a multi-level iteration based on the Landweber-Richardson method for solving moment 
problems. It will turn out that different stopping criteria can be used for stopping the 
Landweber-Richardson method than for the CGNE method. 

The multi-level iteration based on the Landweber-Richardson method is defined as fol- 
lows 

Algorithm 3.1 Let r] > 0, r > 1. 

At level : Without loss of generality we assume that x 5 G X satisfies 

1. ||To(4 )-2/|| >2(l + v )(8+\\P x*-x4), 

2. \\T (x 5 0fi )-y s \\> 2(1 + V )t5. 

If^. does not hold, then Xq is accepted as an approximate solution and the multi- 
level algorithm is terminated. If\^. does not hold, but ^. is satisfied, then we can 
choose a level J such that 

\\Tj(x s 0t0 )-y 5 \\ = ||To« )-/ll >2(l + ri)(6 + \\Pjx,-x.\\) . 

In this situation we simply renumber the levels and set = J. With the new 
setting both and [|. hold, and the "Landweber-Richardson" iteration at level 
is performed 

X k+1,0 = X k,0 ~ ^0 (^0(^,0) ~ y 5 ^) ■ 



s 



The iteration at scale Xq is terminated if for the first time 

11^(4+1,0) -/II <2(l + V )(5+\\P x*-x4) . 
The termination index is denoted by l*(Q, 8). 
From level N to level N + l: If \\T N+1 ( 

x i*(n,5)+i,n) ~ V II — 2(1 + v) T 8> then the multi- 
level algorithm terminates and xf^, NS \ +1 N is an approximate solution. Else we 
define xf, NS \ N =: Xq >n+1 . Since ||P/v+i£* — x* || < ||Pjv^* — x*|| it follows from the 
definition of /* (N, 8) : 

\\Tn+i(4,n+i) ~y S \\ > 2(1 +v) (5 +\\Pn+ix*-x* ||) . 

Therefore, at least one iteration of Landweber-Richardson method at level N+l 
is performed and 

x k+i,N+i = x k,N+i ~ T N+1 (TN + i(x s k N+1 ) — y 5 ^j . 

The iteration at scale X^r+i is terminated if 

\\T N+1 (xi +1>N+1 ) -y s \\< 2(1 +77) (8 + \\P N+1 x* -x4) 

for the first time. The termination index is denoted by l*(N + 1,8). 

In the following, we verify some basic monotonicity properties of the iterates of the 
multi-level algorithm [3.1| . 

Proposition 3.2 Let Assumptions \2.1\ and \2.3j hold. At fixed level N let the iteration 
be stopped according to the generalized discrepancy principle 

(3.1) ||Pjv«jv) - 2/1 > 2(1 + r]){8+\\P N x* -x m \\) for < k < h{N,8) . 
Then 

(3.2) \\xf +1 -P N x*\\ < \\xf -P N x m \\ j = 0,...,h{N,8) . 
Moreover, if (\3. 1\) at level N + l does not become active for xf, N+1 , then 

(3.3) ll x i,iv+i — x *ll 2 — II^jv+i — x *ll 2 — 1 1 „ \\^ nX o,n+i ~ y 5 \\ 2 ■ 
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Proof: From the definition of the multi-level Landweber-Richardson method it follows 



-2(T N (x s k>N - P N x*),T N (4 tN ) - y 5 ) 
(3-4) = \\x s ktN -P N x*f-2\\T N {x{ tN )-y s tf 

+ n (Tn(x{ n ) - y s ) || 2 
+2(T N (P N x*) - y 5 ,T N (x{ N ) - y 5 ) . 



From (|3.1| ) and (|2.2|) it follows 

\\T N (P N x*) -y 5 \\ < 5+ \\T Nt (P N x*) -T Nt (x*) 

(3.5) 

< 5 + \\Pn%* — a:* || • 



For < k < l m (N, 5) it follows from (|J), (gj) and that 



(3-6) +2||T^(P iV x,) -ylHTjv^) -yl 

< \\ x i,N — Pn%*\\ 2 ~ T~7 — \\Tn(xIn) — ^ll 2 



Now ( |3.2| ) can be easily proven with an inductive argument. ( |3.3| ) can be proven with 
similar arguments as above. □ 



From Proposition |37| it follows that ( p.l|) determines a well-defined stopping index for 



the Landweber-Richardson method at level N. 

Corollary 3.3 Let Assumptions \2.1\ and \2.2\ hold. Let the Landweber-Richardson iter- 



ation at fixed level N be stopped according to the generalized discrepancy principle ( \3. 1\ ). 
U 

(3.7) 5+ \\P N x* -x m \\ > 0, 

then 

h(N,6) < oo . 
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Proof: From fl3.1| ) and ( ^.6[ ) it follows 

477(1 + 77) (5 + \\P N x* - xJ\) 2 < 

k=0 



< 
< 

This together with (|3.7|) shows l*(N,5) < 00. □ 

So far we have proven a monotonicity result for the iterates of the Landweber-Richardson 
method at fixed level N. In the following we show that the iterates of the multi- 
level algorithm also satisfy a certain monotonicity when the levels are switched. This 
together with Proposition |3.2| provides monotonicity of the iterates of the multi-level 
algorithm |3.1| . 



\ X U(N,8)+1,N Pn%*\ 
l*(N,6) 



+ 



V 



1 + V fc = 
||-^o,JV Pn%* 
00 . 



£ \\ t n{xI,n) 



.S\\2 



3.1 Convergence and stability of the multi— level algorithm 

Using the monotonicity results of the iterates of the multi-level schemes based on the 
Landweber-Richardson method we verify below that the multi-level scheme is convergent 
and stable. The following lemma will be central for our further considerations. The proof 



of this lemma is analogous to the proof of Lemma 2.5 



Lemma 3.4 Let 6 = 0. If Assumptions 2.1 and 2.i hold, the multi-level algorithm 3.1 



terminates after a finite number of iterations 
or 

there exists a level N such that at each level N < N a finite number of Landweber- 
Richardson iterations are performed and at level N the Landweber-Richardson it- 
eration does not terminate. 

Similarily as in Section ^ (using instead of Hanke's results in the proof of Theorem 
|2.6| some well-known results on convergence properties of the Landweber-Richardson 
method (see e.g. Groetsch [|Gro84| ) it is possible to prove that the multi-level Landweber- 



Richardson algorithm is convergent in the case of exact data. 
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Theorem 3.5 Let Assumptions [O and hold. Then the multi-level Landweber- 
Richardson algorithm converges to a solution x. 



In the case of perturbed data, it follows from Corollary |3.3| that the discrepancy principle 
( |3.1| ) determines a well-defined (finite) stopping index at level N < N, which is denoted 
by /*(iV, 5). In the following we finally terminate the multi-level algorithm if for the first 
time 

(3.8) \\Tn(4,n)-V S \\ <2(l+v)rS. 



From Corollary 3.3 and Assumption |2.2| it follows that this criterion determines a well- 



defined stopping index at level N(S), which is denoted by L*(N(6), S). 

Our next result shows that the multi-level iteration is a stable method. The proof of 
this result is analogous to the proof of Theorem 



Theorem 3.6 Let the assumptions of Theorem \3. f\ hold and assume that Tjy- is contin- 
uously invertible. If y 6 fulfills flg. 3j) , and if the iteration at level N is stopped according 
to the discrepancy principle ( \3.J\) , and the multi-level iteration is terminated by (\3. 8j ), 
then 

X N(S),L t (N(S),S) ~~ * X *i 6^0. 

In practical applications, if the quantity \\P 

Xjf. X ^ 1 1 IS substantially underestimated, it 
may happen that the algorithm does not switch to the next level although the iterations 
stagnate and the approximations do not improve significantly any more. This can be 
avoided for the multi-level iteration based on Landweber-Richardson method by using 
alternatively the following stopping criterion: 

(3.9) \\x% +1 -x%\\> 2(1 + t])t6 for k = 0, l(N, 5) . 

Then, from the definition of the Landweber-Richardson method it follows under the 



Assumption [2J. 



\\(T N x N - y 5 )\\ > \\T* N (T N x N - y s )\\ = \\x% +1 - x%\\ > 2(1 + v )t6 . 

Therefore we have l*(N,8) > l*(N,8). From this it is obvious to see that Proposition 
|3.2| , Corollary |3.3| , Lemma and Theorem [T5] are valid, if we replace in these results /* 
by i*. It is also quite easy to see that (|3J|) determines a well-defined stopping criterion. 



The proof of Theorem |3^ for the multi-level Landweber-Richardson method with the 



stopping criterion |3.9| requires the assumption that is continuously invertible to 



guarantee that the solution of x = y is a solution of T/v*(T^- — y) = and vice versa. 
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4 The multi— level algorithm for the moment prob- 
lem 



In the sequel we show the multi-level algorithms |2.3| and |3.1| can be applied for the 
solution of general moment problems as described in Section [I]. 

It is well-known that a signal x G Xjv can be completely reconstructed from its moments 
(x,gf) if the set {g^}j<zi is a frame for X N . For the readers convenience we quote the 
definition of a frame (see e.g. |You80|| ): 

Definition 4.1 A family {g^}j & j is a called a frame for X^ if there exist constants 
An,Bn, called the frame bounds, such that 

(4.10) A N \\x\\ 2 <J2\(x,gf)\ 2 <B N \\x\\ 2 \fxeX N . 

In this paper (if not stated otherwise) we will always assume that the family {gj}jei is 
a frame for X^. For a given frame {g?}jei we introduce the following linear operator 

(4.11) T N : X N -> , T N x = {(x,gf . 

The operator T N satisfies ||7jv|| < VBn- The adjoint operator is given by 

(4.12) T* N : £ 2 (I) - X N , T* N {^} jeI = . 



We note that if the family {gj}jeJ is an orthonormal basis on Xjf, then Tjy is a unitary 
operator and Tfa = (T/v) -1 . 

As usual we define the frame operator by 

(4.13) S N : X N -> X^, SVx = T*^ = <?f ><?f . 

It is well-known that Sn is positive definite and invertible if {gf}j^j is a frame. Stability 
results for frames in conjunction with moment problems can be found in [Ghr96 . 



Obviously the operator Tn (respectively the scaled operator -jjfT) satisfies Assump- 
tion |2.1| (i) and (iii). A sufficient (but not necessary) condition for the frame analysis 
operators T/v to satisfy Assumption |2.1| (ii) is for instance when the frame elements g^ 



correspond to reproducing kernels, see Section |4.1| below. With the notation introduced 
above this shows the setting in which the generally defined multi-level algorithms of 



13 



Section |2| and Section |3] can be applied to recover x from noisy moment measurements 

For practical applications many delicate problems have to be solved before implemen- 
tation of the multi level algorithms. We need conditions which allow an easy verification 
that the set {gf} is a frame for Ajv, as well as useful estimates of the frame bounds. And 
finally we have to find appropriate estimates for ||Pjvic* — £*|| occurring in the stopping 



criteria (p.4|) and (|3.1| ) 



In Section 5 we summarize recent results from the literature which allow to estimate 
frame bounds efficiently 

4.1 Signal recovery in reproducing kernel Hilbert spaces 

A special case of the moment problem occurs in reproducing kernel Hilbert spaces. Let 
A be a Hilbert space of real or complex valued functions defined on a subset I of R with 
inner product (•, •). We recall that a Hilbert space A is called a reproducing kernel Hilbert 
space if for every element t G I the "point evaluation" x \— > x(t) on A is bounded. From 
the Riesz representation theorem it follows that there exists a unique K t G A, which 
satisfies 

(4.14) (x, K t ) = x(t) for every x G A 
The reproducing kernel K on / x / is defined by 

(4.15) K(t,u):=(K t ,K u ) = K t (u) . 

In this section for application of the multi-level algorithm of Section Q we use a family 
of reproducing kernel Hilbert spaces Ajv- Note that every closed subspace of a RKHS 
is itself a RKHS. Thus a reproducing kernel Hilbert space always contains a family of 
nested reproducing kernel subspaces. 

If A^r C A is a RKHS, then the orthogonal projector of A onto Ajy can be charac- 
terized by (see e.g. [|Ist87| , Chapter 11]) 



(4.16) {P Xs x)(t) = {x,K»). 

The following lemma is an immediate consequence of the above characterization of the 
orthogonal projector P N 

Lemma 4.2 Let xm G Xm, where M > N. Moreover, let Tjy be defined by 

Tn(-) := (;9?), 
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where = K N (-,j). Then 

(4.17) T N (P N x M ) = T N x M . 



Using Lemma |4.2| we are able to verify the general assumptions in the convergence 
and stability results of the multi-level algorithm. Moreover, we discuss details of the 
numerical implementation, especially implementation of the stopping criteria. 



For reproducing kernel Hilbert space the stopping criterion (|3.1| ) is equivalent to 

\J24A^)-y 5 \ < 2(1+77) (6 + \\p n x*-x4) . 

There are several important examples of reproducing kernel Hilbert spaces which 
deserve to be mentioned explicitly. 

(i) Let x G Bn where Bn is the space of band-limited functions, i.e. 

/oo 
\x(t)\ 2 dt < oo and supp£ C [-N, N] , 
-oo 

where x denotes the Fourier transform of x. The space Bn is a RKHS with reproducing 
kernel 

(4.19) K N (t,u) 



rN A sin N(t - u) 



N(t - u) 



This space is of paramount importance for digital signal processing, see e.g. ||Hig96|| . We 
will deal with Bn in more detail in Section [5] in conjunction with signal recovery and 
bandwidth estimation from nonuniformly spaced noisy samples. 

(ii) The RKHS Vn consists of all trigonometric polynomials p of degree N given by 

N 

(4.20) p(t)= £ a n e 2mnt . 

n=-N 

The reproducing kernel is K(t, u) = DN{t — u) where the Dirichlet kernel is 

sin(2iV+ l)7rt 



(4.21) D N (t) 



sin Ttt 



It follows immediately from the fundamental theorem of algebra that the set {gj}j£i, 
where gj(t) = D^it — tj) constitutes a frame for Vn whenever |/| > 2N ', where the sam- 
pling points tj may be randomly spaced. Grochenig has derived estimates of the frame 
bounds for the weighted frame {wjgj}j e i, if the maximal gap between two successive 
sampling points satisfies a certain criterion, see ||Gro93 . 



(iii) Sobolev spaces H r for r > | are well known to possess reproducing kernels. Note 
however that the inner product on H r is not the restriction of that of L 2 (IR). Hence the 
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formulae in the previous sections must be reinterpreted in terms of the H r inner product. 
See |[N W91|| for a detailed discussion of Sobolev spaces and reproducing kernels. 



Wavelet subspaces that constitute a so-called multiresolution analysis (MRA) ||Mal89 



Dau92] are further important examples for reproducing kernel Hilbert spaces [[XN94 1 



where moment problems arise naturally when one wants to reconstruct a function from 
its wavelet coefficients. Recently the concept of a MRA has been generalized, leading to 
a so-called frame multiresolution analysis |BL97| . More examples for reproducing kernel 



Hilbert spaces can be found in |[N W91| , [lst87|| 



5 Reconstructing a signal of unknown bandwidth 
from irregularly spaced noisy samples 

In many practical applications, for instance in geophysics, spectroscopy, medical imaging, 
the signal cannot be sampled at regularly spaced points. Thus one is confronted with the 
problem of reconstructing an irregularly sampled signal. Often the signal can be assumed 
to be band-limited (or at least essentially band-limited) , but the actual bandwidth is not 
known. 

The problem of reconstructing a band-limited signal x from regularly and irregularly 
samples {x(tj)}j & i has attracted many mathematicians and engineers. Many theoretical 
results and efficient algorithms have been derived in the last years, see e.g. ||Zay93| , 



|FG94j , |Sam95| , |Hig96|| and the references cited therein. In ||Chr96|| results on frame 



perturbations are applied to the irregular sampling problem to derive error estimates 
for the approximation. Most of these results are based on the assumption that the 
bandwidth of the signal to be recovered is known a priori. However in many applications 
this assumption is not justified. This is the motivation to consider the problem of 
reconstructing a band-limited signal whose actual band-width is not known, given only 
a set of nonuniformly located and noisy samples. 

In the sequel we show that the multi-level algorithm developed in Section [] in com- 
bination with the results of Section [| provides an efficient solution to this problem. To 
do this, we verify the general assumptions in the convergence and stability results of the 
multi-level algorithm. 

In the approach suggested below, we use the adaptive weights method developed by 



Feichtinger and Grochenig [FG94]. The adaptive weights method is based on the fact 
that a signal x G Bn can be reconstructed from its nonuniformly spaced samples, if the 
maximal gap between two successive samples does not exceed the Nyquist interval tt/N. 
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The precise formulation of this results is summarized in the next theorem 



Theorem 5.1 Let {tj} ^£ be a strictly increasing sequence of real numbers (called the 
sampling set), which satisfies 

(5.22) n/N > 7 := sup(t j+1 - tj) . 

"N 1 



Then {^JuTjK^} is a frame for with frame bounds 

(5.23) A= (1-— ) 2 andB = (1 + 2^) 2 . 

7T 7T 

where the weights Wj are given by 

(5.24) w j =(t j+1 -t j - 1 )/2. 



It is a consequence of Theorem [5J] that a signal x £ Bjy can be recovered from its 
sampled data {x(tj)} by the Landweber-Richardson method, where the operators T N 
and are characterized by 

(5.25) T N x = {(x, K^)^/w]} jeZ , T^{ Cj } ja = £ ^/w]CjK^ . 

In the remaining of this section, we assume that for given sampling set {^} Jg ^ the family 
{y/WjK^}^^ satisfies (|5.22|) and thus forms a frame for B^. 

For an efficient implementation of the multi-level algorithm practical relevant esti- 
mates for ||Pjv^* required. In this example Pn denotes the orthogonal projec- 
tor of L 2 onto Bn. Assume that the noisy samples xl(tj) := x*(tj) + 5j are given with 
SjeZ I I = ^ an d x * e &m f° r fixed M. It holds 

llx* - P/v£*|| 2 = ||x*|| 2 - HP/vx*!! 2 = / \x*(u)\ 2 du> . 

J\w\>N 

for N < M, where we have used Plancherel's equality. Here x denotes the Fourier 
transform of x. In the sequel we show how to estimate — P^x^ 2 recursively. 
It follows from the definition of frames that 

(5.26) \M t j)\ 2w j < \\ x *\\ 2 < -i- \ x *(tj)\ 2w i 

with equality if the sampling points are regularly spaced. In this case it is easy to see that 
{y/WjK^}^^ is a tight frame for Bn for N < M. Note that |5.26| is only a worst-case 
estimate. Extensive numerical experiments have shown that for practical applications 

(5.27) EWfe)lH- W 2 «« 
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as long as T is not very ill-conditioned. 

If no a priori information about ||P/v£* — P/v_i£*|| 2 is available, a reasonable estimate 
for — P/vx*|| 2 is given by — P/v-iX* || 2 - Let us start with the first level, i.e. N — 1. 
In this case we have \\x* — Pia:*|| 2 < — Poa;*|| 2 = II s * II 2 - By virtue of (|5.27| ) we can 
use g ^ \ x l(~tj)\ 2w j as estimate for \\x* — Pia;*|| 2 . Now we run algorithm [2]^ or |3J] at 



the first level until the corresponding stopping criterion (|2.4j ) or (|3.1|) is satisfied. We 
obtain an approximation denoted by x 5 kt 1 and use 1 || 2 as approximation to ||P 2 a; !| ,|| 2 
to estimate ||x* — P2£*|| 2 in the next level. We proceed inductively by using \\x 5 kt ^H 2 to 
estimate 

(5.28) \\x* - Pv +1 x,|| 2 « \4{tj)\ 2 Wj - lK,ivl| 2 

used in the stopping criterion at level N + 1 . 

Improved estimates for ||x* —Pjyx* || 2 can be obtained by trying to estimate ||P/v_i2* — 
P/v^*|| 2 = S\u\=n-i |^*(^)| 2 ^ using a priori information about the decay of the frequen- 
cies of x(u). This can be done either by a model-based approach or by a statistical 
analysis of a given class of signals. For instance potential fields, such as the gravity or 
the magnetic field, have exponentially decaying frequencies IjRSaifl . It is well-known that 
the Fourier coefficients of many bio-medical signals and large classes of images decay like 
l/u a where a is a constant depending on the class of signals under consideration. 



6 Experimental results 

We demonstrate the performance of the proposed algorithms for the reconstruction of a 
band-limited signal from its nonuniformly spaced noisy samples. For the implementation 
of the algorithm we use the discrete model presented in |FGS95|[ . 



Experiment 1: In the first example we consider a signal from spectroscopy. The original 
signal x* of bandwidth M = 30 consists of 1024 regularly spaced and nearly noise-free 
samples. We have added zero-mean white noise with a signal-to-noise ratio of 12% 
and have sampled this noisy signal at 107 nonuniformly spaced sampling points. A 
simple measure for the irregularity of the distribution of the sampling points is the 
standard deviation a of the differences of the sampling points. Let Aj = tj + \ — tj and 
s = Y7j=i Aj/r then a is given by 



a 



/ s; = iiA 3 -- a | 

r — 1 



2 
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Obviously a = for regularly sampled signals. In this experiment the maximal gap 7 
in the sampling set is 33 and a = 6.3. The setup of this experiment simulates a typical 
situation in spectroscopy. 

Since we have to reconstruct x* from noisy samples, we cannot expect to recover 2* 
completely. 

In our experiments we have not used any a priori information about the decay of the 
frequencies of the signal, but have estimated \\x* — Pnx*\\ by the simple recursive scheme 
as described in the previous section. A comparison of \\x* — Pjyx*\\ 2 and J2j \ x l(tj)\ 2w j ~ 
\\ x t* n\\ 2 i s demonstrated in Figure |6TT| . Algorithm [D] terminates at level iV* = 19, the 



reconstruction is shown in Figure |T2] (solid line). The dashed line represents the original 
signal, and the * are the noisy samples. The normalized error x^. tv 1 1 / 1 1 1 1 ^ 

0.17. 

Recall that we have had to use a different stopping criterion than for Algorithm |3.1| 
to prove convergence of the CG-multi-level algorithm. The CG-stopping criterion (|2.4|) 
forces the algorithm to terminate earlier, thus giving a weaker approximation, see Fig- 
ure |6.3j . Although we cannot justify it theoretically, we always got very good recon- 
structions in our numerical experiments by using the same stopping criterion as for the 
Landweber- Richardson multi- level algorithm, compare Figure ET1. 



Note that our numerical experiments indicate that for most practical applications 
one can drop the factor 2 in the discrepancy principle (|3.8f) . In this case Algorithm |3J] 
terminates at level 23 and the approximation error is 0.09, see also Figure IO. 



Experiment 2: We use the same signal and the same noise level as in Experiment 1, 
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gure 6.2: Approximation using Algorithm |3.1| , terminated at level 19 
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jure 6.3: Approximation using Algorithm |2.3| , terminating at level 13 



20 







1 




sk Noisy samples 








— Original signal 








Approximation 

















100 200 300 400 500 600 700 800 900 1000 



Figure 6.4: Approximation using Algorithm [2.3|, using stopping criterion [3.1|, terminating 
at level 19. 
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Figure 6.5: Approximation using Algorithm |3.1] , omitting the factor 2 in the stopping 
criterion |3.1|, terminating at level 23. 



21 



x Noisy samples 
— Original signal 
Approximation 




100 200 300 400 500 600 700 800 900 1000 



Figure 6.6: Approximation using Algorithm |3J] providing an approximation error \\x* 
x\\ = 0.19. 



x fc»,AT, 



but only 89 sampling points that are in addition more irregularly distributed with a 
standard deviation o = 12.2. Thus the operators T/v have a larger condition number 
than those arising in Experiment 1. This experiment demonstrates that a proper choice 
of the bandwidth can serve as a regularization procedure. Algorithm |3.1| terminates at 
level A* = 18, the final error between the original signal and the approximation is 0.19. 



The reconstruction is shown in Figure |6.6j . It is remarkable that this approximation is 
better than the approximation we obtain when we apply the standard CGNE method E7T 
directly using the a priori information about the "correct level" M = 30. In this case 
the approximation error is 0.28, thus substantially larger. The approximation is shown 
in Figure $7f[ 



Experiment 3: We consider the reconstruction of a signal representing a dynamical 
response of a mechanical system. The signal has been sampled at 300 irregularly spaced 
points with a = 40.7, the original signal itself is not known. We run Algorithm |3.1| 
without using any other a priori information about the signal than the samples and the 
measurement error which is known from experience to be approximately 10%. 

Algorithm terminates at level 42, the approximation is displayed in Figure o^E as 
solid line. In the chosen discrete model IIFGS95I, N = 42 means that the reconstruction 



is a trigonometric polynomial of degree 42. Since the number of sampling values is 300, 
we have about 3.5 times as many samples as we would need theoretically to determine 
a reasonable approximation. Thus we should obtain about the same bandwidth and 
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Figure 6.7: Approximation using standard CGNE directly for level M = 30, providing 
an approximation error \\x* — x 1 ^ II /II x II = 0.28. 

approximation by using only every second sampling value, since we still have enough 
samples and the corresponding operators Tjy are still well-conditioned. It turns out that 
the algorithm terminates again at level 42, giving an approximation close the first one. 
The normalized error between both approximations is 11%. Having in mind that we 
have used only 50% of the given samples and that the measurement error is about 10%, 
this result demonstrates the robustness of the proposed algorithms. The approximation 
for using only half of the samples is displayed in Figure |6.8| as dotted line. 
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Figure 6.8: Two reconstructions of dynamical impulse response using Algorithm 3.1 
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